class: inverse,left, middle background-image: url(data:image/png;base64,#background.png) background-size: cover <img src="data:image/png;base64,#LOGO_DIPLOMADO.png" width="500px"/> ##Módulo 3: Percepción remota, análisis masivo y GEE ### Introducción a las Series de tiempo Matías Olea <br> <a href="https://orcid.org/0000-0003-0194-7784"> ORCID </a><br> matias.olea@pucv.cl</a><br> .large[<b><a href="https://www.pucv.cl/uuaa/site/edic/base/port/labgrs.html">LabGRS</a> | Noviembre 2023</b>] <br> --- class: center,middle background-image: url(data:image/png;base64,#labgrs_logo.png) background-size: 35% --- ## Contenidos .pull-left[ 1) Que es una serie de tiempo. 2) Extracción de series de tiempo (vector) desde raster temporales. 3) Rellenado de datos. 4) Remuestreo temporal. 6) Componentes de una serie y descomposición. 7) Análisis de cambios en tendencias (ciclos). 8) Anomalías respecto a desviaciones de la media (Z-score). ] .pull-right[ <img src="data:image/png;base64,#https://raw.githubusercontent.com/allisonhorst/stats-illustrations/main/rstats-artwork/r_rollercoaster.png" width="650px"/> ] --- ### Series de tiempo Una **serie de tiempo** corresponde a un conjunto de observaciones de variables cuantitativas en el tiempo. Estas variables pueden o no tener distintos comportamiento en el tiempo.
--- ### Series de tiempo Las series las podemos encontrar como estacionarias y no estacionarias .pull-left[ **Estacionaria**. Es aquella cuya media y varianza se mantienen en el tiempo, es decir, son constantes, por lo que facilita más su modelamiento.
] .pull-right[ **No Estacionaria**. Es aquella cuya media y varianza cambian en el tiempo, es decir, no son constantes, por lo que dificulta más su modelamiento.
] --- ### Series de tiempo Las series temporales poseen 4 componentes (3 para algunos autores) los cuales son en simples palabras: -- **Tendencia** (T - trend): corresponde al comportamiento de la media y la varianza en el tiempo. -- **Ciclos** (C - cycle): oscilaciones respecto a la tendencia en largas o medias escalas temporales. Como se desprende de la tendencia, algunos autores no la consideran como un componente propiamente tal. -- **Estacionalidad** (S - seasonality): comportamientos determinados en periodos fijos. Corresponden a cambios dentro de unidades de tiempo de manera regular (como cambios intra-anuales). -- **Aleatorio** (R - random): comportamiento sin un patrón determinado. Es un comportamiento practicamente impredecible. Lo podemos encontrar en la literatura también como "erratico". --- ### Series de tiempo <center><img src="data:image/png;base64,#raster_ts.jpg" width="1000px"/></center> --- ### Series de tiempo #### Lectura de datos ```r library(terra) library(tidyverse) ``` ```r setwd("C:/TuDirectorio/") # Leemos fechas desde 2001 y /10000 para pasar de numeros enteros a decimales myd_EVI <- rast("MYD13Q1_EVI_ts491.tif")/10000 # Dividimos por 10000 por escalado del producto tab_dates <- read_csv("MYD13Q1_dates_ts491.csv") # leemos nuestra tabla con fechas ``` ```r nlyr(myd_EVI) # número de bandas (1 banda = 1 fecha) ``` ``` ## [1] 491 ``` ```r nrow(tab_dates) # número de fechas ``` ``` ## [1] 491 ``` --- ### Series de tiempo #### Creación serie de tiempo a partir de 1 pixel ```r xy_plant <- cbind(209482, 6085390) %>% as.matrix() # Coord. Plantación con incendio 2017 xy_nativo <- cbind(209362, 6041365) %>% as.matrix() # Coord. Bosque nativo deciduo ``` ```r serie_plant <- terra::extract(myd_EVI, xy_plant) %>% as.numeric() # Extraemos serie 1 pix serie_nativo <- terra::extract(myd_EVI, xy_nativo) %>% as.numeric() ``` ```r tab_pixel <- tibble(dates = rep(tab_dates$date, 2), EVI = c(serie_plant, serie_nativo), Serie = c(rep("Plantacion Mixto", 491), rep("Nativo Deciduo", 491))) ``` ```r write.csv(tab_pixel, "Tabla_muestreo_pixel.csv", row.names = F) ``` .footnote[ <span style="font-size:9pt"> Como vimos en clases anteriores, en lugar de usar una matriz con las coordenadas XY de un pixel, podemos utilizar un polígono incorporando el estadístico que queremos usar y si queremos incluir los pixeles que toncan los bordes del políno con el argumento touches </span> ] --- ### Series de tiempo #### Creación serie de tiempo a partir de 1 pixel <img src="data:image/png;base64,#DIPGEOPR_03_9_files/figure-html/unnamed-chunk-12-1.png" width="100%" /> --- ### Series de tiempo #### Rellenado de la serie Existen muchos paquetes que nos permiten realizar rellenos o suavizado de nuestros datos como **gapfill**, **oce**, **imputeTS**, **TSstudio**, o **phenofit**, entre otros. Los más frecuentes para análisis de datos de IV (Índice de Vegetación) proveniente de datos satelitales es **phenopix** o el software **TIMESAT**. Otra opción es la creación de Composiciones (Composites), que realizan agregaciones temporales en ciertos intervalos de tiempo para asegurar una mayor cantidad de datos disponible. Un claro de ejemplo de esto es el producto MOD13Q1 de Modis que realiza composites de 16-días. Para este ejercicio ocuparemos el paquete imputeTS. ```r library("imputeTS") ``` --- ### Series de tiempo #### Rellenado de datos ```r tabla_plantacion <- tibble(dates = tab_dates$date, EVI_raw = serie_plant) ``` Para las funciones de imputeTS, necesitamos un objeto clase Time-Series: ```r ts_plant <- ts(serie_plant, tab_dates$date) ``` ```r tabla_plantacion$Lin <- na_interpolation(x=ts_plant, # nuestro objeto TS option = "linear", # metodo de rellenado maxgap = Inf) %>% # n de NA consecutivos para rellenar as.numeric() # resultado a vector numerico para incluirlo en la tabla ``` --- ### Series de tiempo #### Remuestreo temporal El remuestreo temporal consiste en realizar **agregaciones** o **agrupaciones** temporales, según una condición de tiempo (mensual, anual, semestral, trimestral, etc.). Comunmente, estos ejercicios se realizan para: - Tener una **periodicidad** o **intérvalos de tiempo regulares** en nuestros datos. Ej: Caso Landsat (16-días 1 satelite / 8-días 2 satelites) - **Disminuir las brechas** por pérdida de datos causados por la presencia de nubes. Ej: Composites de MODIS. De datos diarios a cada 16-dias. - Homologar fuentes de información con **diferentes resoluciones temporales**. Ej: Datos de dendrocronológicos (1-año) con datos satelitales o climáticos. --- class: center,middle background-image: url(data:image/png;base64,#tiempo.png) background-size: 90% background-color: black --- ### Series de tiempo #### Remuestreo temporal Para este ejemplo, utilizaremos nuestra serie de tiempo EVI de MODIS-Aqua de 16-días (23 datos al año) para remuestrearla temporalmene en datos mensuales. **Opción 1: usando tidyverse** ```r tabla_plantacion$month <- month(tabla_plantacion$dates) # creamos una nueva columna con los meses tabla_plantacion$year <- year(tabla_plantacion$dates) # creamos una nueva columna con los años tabla_mensual <- tabla_plantacion %>% group_by(year, month) %>% # agrupamos por año y mes summarise(EVI_mon = mean(Lin, na.rm = T)) # agrupamos usando la media tabla_mensual$dates <- paste(tabla_mensual$year, tabla_mensual$month,"15", sep = "-") %>% as.Date() # construimos vector de fechas para las agrupaciones ``` --- ### Series de tiempo #### Remuestreo temporal <img src="data:image/png;base64,#DIPGEOPR_03_9_files/figure-html/unnamed-chunk-18-1.png" width="100%" /> --- ### Series de tiempo #### Remuestreo temporal **Opción 2: usando xts** ```r library(xts) # la versión para raster es "rts" xts_plantacion <- xts(tabla_plantacion$Lin, # vector numerico con la variable tabla_plantacion$dates) # vector de fechas EVI_mensual <- apply.monthly(xts_plantacion, FUN = mean) # agregación mensual ``` --- ### Series de tiempo #### Remuestreo temporal <img src="data:image/png;base64,#DIPGEOPR_03_9_files/figure-html/unnamed-chunk-20-1.png" width="100%" /> --- ### Series de tiempo #### Test de estacionariedad ```r library("tseries") ``` **Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test** **H<sub>0</sub>**: La serie de tiempo tiene una tendencia estacionaria. **H<sub>A</sub>**: La serie de tiempo tiene una tendencia no estacionaria. --- ### Series de tiempo #### Test de estacionariedad ```r kpss.test(tabla_plantacion$Lin, null = "Trend", lshort = F) # Plantacion Mixto (rellenada) ``` ``` ## ## KPSS Test for Trend Stationarity ## ## data: tabla_plantacion$Lin ## KPSS Trend = 0.30546, Truncation lag parameter = 17, p-value = 0.01 ``` Se rechaza H<sub>0</sub> por lo tanto la serie es no estacionaria. <img src="data:image/png;base64,#DIPGEOPR_03_9_files/figure-html/unnamed-chunk-23-1.png" width="100%" /> --- ### Series de tiempo #### Test de estacionariedad ```r kpss.test(na.omit(serie_nativo), null = "Trend", lshort = F) # Nativo (con NA) ``` ``` ## ## KPSS Test for Trend Stationarity ## ## data: na.omit(serie_nativo) ## KPSS Trend = 0.10429, Truncation lag parameter = 17, p-value = 0.1 ``` No se rechaza H<sub>0</sub> por lo tanto la serie es estacionaria. <img src="data:image/png;base64,#DIPGEOPR_03_9_files/figure-html/unnamed-chunk-25-1.png" width="100%" /> --- ### Series de tiempo #### Descomposición de la serie ```r ts_mensual <- ts(tabla_mensual$EVI_mon, # vector con la variable start = c(2002, 7), # año y mes de inicio de la serie end = c(2023, 10), # año y mes de termino de la serie frequency = 12) # frecuencia de datos al año dec_plant <- decompose(ts_mensual) # descomposición ``` --- ### Series de tiempo <img src="data:image/png;base64,#DIPGEOPR_03_9_files/figure-html/unnamed-chunk-28-1.png" width="100%" /> --- ### Series de tiempo #### Análisis de cambios en tendencias (ciclos) ```r library(DBEST) ciclos <- DBEST(data = ts_mensual, # objeto clase Time-Series data.type = "cyclical", # datos con o sin componente estacional seasonality = 12, # cantidad de datos por ciclo anual algorithm = "change detection", # tipo de detección breakpoints.no = 3, # numero de breaks solicitados para encontrar first.level.shift = 0.05, # umbral principal de cambio second.level.shift = 0.1, # umbral secundario de cambio distance.threshold = "default", # busqueda de umbral si está por defecto duration = 24, # ventana de tiempo para busqueda alpha = 0.05) # p-value ``` --- ### Series de tiempo <img src="data:image/png;base64,#DIPGEOPR_03_9_files/figure-html/unnamed-chunk-31-1.png" width="100%" /> --- ### Series de tiempo ```r ciclos$Start # [primer, segundo y tercer más grande] ``` ``` ## [1] 174 187 1 ``` ```r ciclos$End # [primer, segundo y tercer más grande] ``` ``` ## [1] 175 202 28 ``` ```r tabla_mensual$dates[ciclos$Start[1]] ``` ``` ## [1] "2016-12-15" ``` ```r tabla_mensual$dates[ciclos$End[3]] ``` ``` ## [1] "2004-10-15" ``` --- ### Series de tiempo #### Anomalías respecto a desviaciones de la media (Z-score) usando periodo completo ```r tabla_mensual$Zscore <- scale(tabla_mensual$EVI_mon, scale = T) ``` <img src="data:image/png;base64,#DIPGEOPR_03_9_files/figure-html/unnamed-chunk-34-1.png" width="100%" /> --- ### Series de tiempo #### Anomalías respecto a desviaciones de la media (Z-score) usando ciclo más estable ```r tabla_mensual$Zscore_2 <- (tabla_mensual$EVI_mon - mean(tabla_mensual$EVI_mon[28:174]))/sd(tabla_mensual$EVI_mon[28:174]) ``` <img src="data:image/png;base64,#DIPGEOPR_03_9_files/figure-html/unnamed-chunk-36-1.png" width="100%" /> --- ### Bibliografía (2012) E. B. Brooks, V. A. Thomas, R. H. Wynne and J. W. Coulston, "Fitting the Multitemporal Curve: A Fourier Series Approach to the Missing Data Problem in Remote Sensing Analysis," in IEEE Transactions on Geoscience and Remote Sensing, vol. 50, no. 9, pp. 3340-3353. doi: 10.1109/TGRS.2012.2183137. (2004) Jin Chen, Per. Jönsson, Masayuki Tamura, Zhihui Gu, Bunkei Matsushita, Lars Eklundh, A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter, Remote Sensing of Environment, Volume 91, Issues 3–4, (2022) Kong, D., McVicar, T. R., Xiao, M., Zhang, Y., Peña-Arancibia, J. L., Filippa, G., Xie, Y., Gu, X. phenofit: An R package for extracting vegetation phenology from time series remote sensing. Methods in Ecology and Evolution, 13, 1508– 1527. https://doi-org.pucv.idm.oclc.org/10.1111/2041-210X.13870 (2015) G. Yang, H. Shen, L. Zhang, Z. He and X. Li, "A Moving Weighted Harmonic Analysis Method for Reconstructing High-Quality SPOT VEGETATION NDVI Time-Series Data," in IEEE Transactions on Geoscience and Remote Sensing, vol. 53, no. 11, pp. 6008-6021. doi: 10.1109/TGRS.2015.2431315. --- class: inverse middle 